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Abstract. The focusing of acoustic waves is used to study nucleation phenomena in liquids. At large am- 
plitude, non-linear effects are important so that the magnitude of pressure or density oscillations is difficult 
to predict. We present a calculation of these oscillations in a spherical geometry. We show that the main 
source of non-linearities is the shape of the equation of state of the liquid, enhanced by the spherical geom- 
etry. We also show that the formation of shocks cannot be ignored beyond a certain oscillation amplitude. 
The shock length is estimated by an analytic calculation based on the characteristics method. In our nu- 
merical simulations, we have treated the shocks with a WENO scheme. We obtain a very good agreement 
with experimental measurements which were recently performed in liquid helium. The comparison between 
numerical and experimental results allows in particular to calibrate the vibration of the ceramics used to 
produce the wave, as a function of the applied voltage. 
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1 Introduction 
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Recent experiments have shown that acoustic waves can 
be used to study the nucleation of phase transitions far 
from equilibrium under very clean conditions |[|^] . Thanks 
to hemispherical piezo-electric transducers, we have fo- 
cused IMHz acoustic waves in liquid helium and produced 
large pressure and density oscillations. These waves are 
quasi-spherical and, at the acoustic focus (the center), 
their amplitude can be very large. We used an optical 
method to detect the nucleation of bubbles by the neg- 
ative swings of the waves [Q,^. This nucleation occurs 
beyond a certain threshold in the sound amplitude which 
needs to be determined as accurately as possible, in order 
to compare with independent theoretical predictions. We 
later obtained evidence for the nucleation of crystals by 
the positive swings and had the same need. 

In the absence of non-linear effects, the measurement 
of the nucleation threshold would be simple to do. For 
example, one could study the nucleation as a function of 
the static pressure in the experimental cell, and then use 
a linear extrapolation However, we expect non-linear 
effects to occur, especially in cavitation studies. Indeed, 
the homogeneous nucleation of bubbles occurs near the 
"spinodal limit" where the compressibility diverges and 
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the sound velocity vanishes. When an acoustic wave is 
produced in a fluid, with an amplitude such that during 
the negative swings the sound velocity approaches zero, it 
is clear that the wave must be highly distorted. Non-linear 
effects have been already noticed by several authors [^,||, 
§■ 

We thus try to calculate the non-linear focusing of 
the acoustic waves. We start with the spherical geome- 
try, because in a first approximation, everything depends 
only on the radial distance r from the center. As we shall 
see (section H), this calculation still appears difficult be- 
cause the focusing of acoustic waves leads to the forma- 
tion of shocks at all amplitudes in a spherical geometry, 
and their treatment is not quite straightforward. We first 
obtain this result and the associated shock length from 
an analytic calcu lati on which uses the methods of charac- 
teristics (section 2^). Our calculation extends the former 
work of Nemirovskii |^ to the spherical case, except that 
we neglect the coupling with heat modes. It is done in 
the spirit of Greenspan and Nadim's work j7j, though in 
our case it is slightly more tricky due to the shape of the 
equation of state. We make it quantitative by using the 
equation of state of liquid helium |^ which is well estab- 
lished. For weak oscillation amplitudes of the transducer, 
shocks can be ignored and the pressure calculated at the 
focal point by simulating the Euler equations using a fi- 
nite difference method (section 2.3). Indeed, shocks form 
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with an infinitesimal amplitude at a distance r from the 
center which is much less than our mesh size, so that one 
can neglect them. At larger amplitude, shocks cannot be 
ignored. In order to treat the shocks, we have adapted to 
the case of helium a code devoted to shock simulation, 
based on a WENO scheme (section |2^ ). In the end, we 
obtain the amplitude of the density oscillation at the focus 
as a function of an important parameter, the amplitude 
of the displacement of the transducer surface where waves 
are generated. This transducer is a piezo-electric ceramic. 
We choose IMHz for the frequency of the waves, in or- 
der to compare with the experiments. In parallel, we have 
built an experiment to measure the focusing in a quasi- 
spherical geometry. As explained in section |3[ the results 
of this experiment allow a precise comparison with our 
theoretical and numerical work We find that the 

shape of the acoustic wave is indeed distorted at high am- 
plitude and very well described by our calculations, thus 
validating our theoretical method. The final comparison 
with our calculation allows us to calibrate the efficiency 
of the ceramics. As described in our conclusion, this work 
should now be extended to different geometries. One of 
them is the hemispherical geometry where, according to 
other experimental results (2|, non- linear effects are ap- 
parently less important, an observation which needs to be 
understood and compared with future calculations. 



2 Theory 

2.1 Description of the model 

Throughout this paper, we consider a spherical geome- 
try. We take it as one-dimensional since the pressure and 
density fields only depend on the radial distance r. We 
neglected dissipation since our main goal was to compare 
with experiments in superfluid helium 4 which has zero 
viscosity and where the attenuation of sound vanishes in 
the low temperature limit. 

In the case of liquid helium 4 at zero temperature, 
the equation of state has been obtained by three differ- 
ent methods (sound velocity extrapolations, density func- 
tional calculations, and Monte Carlo simulations) with 
similar results. Maris uses the simple form to relate 
the pressure (P) to the density {p): 



with 



Psp = -9.6435 bar 

-3 



Psp = 94.18 kg m 
h 



Ao CsqT = 0.238 mm for 1 MHz waves whose period is 

r = 1// = 1/is. 

We also considered helium 3, a lighter liquid which is 
not superfluid except at very low temperature - i.e. below 
the achievable temperature in our experiment. The same 
form is used for the equation of state, now with 

Psp = -3.1534 bar 
Psp — 53.50 kg m^"^ 
h = 19.262 m-^s^^kg^^ 

The value of Pgp is less negative, which means that the 
inner cohesion of liquid helium 3 is weaker than for he- 
lium 4. At P = 0, the density of liquid helium 3 is po = 
81.916 kgm^"^. The sound speed is then Cgo = 182. 5ms^^ 
and the wavelength Aq — CsqT — 0.182 mm. 

At the temperatures considered, the viscosity is very 
weak and can be neglected though, in case of helium 3, 
neglecting dissipation is an approximation which would 
need to be better justified. Therefore, the numerical ap- 
proximation considers the Euler equations. In order to use 
a dimensionless form of the equations, we have chosen as 
a time scale the period of the wave T, as a length scale the 
wavelength at zero pressure Aq = CsqT and as a density 
scale the spinodal density psp- If we now consider p as the 
dimensionless density and u as the dimensionless velocity 
of the fluid, the Euler equations are written as follows: 



dtp + udrP + pdrU 
dtu + udrU 



-Opu 



drP, 



(2) 



or, by using the conservative variables p and j = pu, 

~9j 



dtp + drj = 

dtJ+dr[f/p + P{p)] = 



-Of 

pr 



(3) 



where 9 is respectively 0, 1, or 2 in planar, axisymetric 
cylindrical, or spherical geometry. The equation of state, 
governing the dimensionless pressure variations, reads: 



P (^2 

PspCgO J 



(4) 



14.030 m^s^^kg"^ 



with Co = 1.848 for hehum 4 and Co = 1.883 for helium 3. 
Note that, in order to have simple notations, we use the 
same names for physical and reduced variables. In gen- 
eral, calculations will be performed with reduced variables, 
while numerical results and experimental parameters will 
be given as physical quantities. 

Boundary conditions are imposed at the center (r — 0) 



Psp is the spinodal limit where the compressibility di- 
verges, the sound velocity vanishes and the liquid be- 
comes totally unstable against the formation of the va- 
por. At P = 0, the density is po = 145.13 kg m"-''. The 
sound speed is then Cso = 238.3 ms^^ and the wavelength 



u = 0, 



(5) 



and on the transducer surface {r = Lq) 

u{t) = —lli Axq sin(ajt) [1 — exp(— 1/1.5)] . (6) 
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Indeed, the motion of the transducer surface is much smaller 
than a mesh size in all our simulations, and it is sufficient 
to impose an oscillating velocity on the ceramics. The ex- 
ponential term represents the response time of the trans- 
ducer. In the experiment the response time is rather equal 
to 8 fis. Here we took it shorter to have a more rapid con- 
vergence of the calculation to the steady regime. This has 
no effect on the final result, as our simulations will always 
be used in the stationary regime. 

In the whole paper, all the calculations are based on 
Euler equations, for which shock waves do occur. We are 
aware that, in superfluids as helium 4, there is no shock 
wave strictly speaking. Actually, when a steep gradient ap- 
pears, it is regularized by dispersion instead of dissipation. 
This means that modes propagate with a different velocity 
depending on their frequency, and thus highest frequency 
modes are expelled from the steep region. However, we 
assume that, as shocks form only in the focal region, and 
during a limited time, not too much momentum is lost 
locally due to dispersion, and that in a first approxima- 
tion, a viscous or dispersive regularization of the shocks is 
equivalent (note that any numerical scheme able to han- 
dle shocks will always introduce either a dissipation or a 
dispersion term in order to make shocks regular, even if 
Euler equations do not contain any viscous term). Besides, 
as we shall see (section 2.3.2), we do not need to describe 
the shock structure exactly, as we are above all interested 
in the relaxation part of the wave. We have found from 
our simulations that the negative pressure swing only de- 
pends very weakly on the numerical viscous regularization 
we use - and thus on the shock amplitude near the focal 
point. Still, it would of course be interesting to be able to 
quantify further the effects of dispersion, and this could 
be the object of a further work. 

Throughout the paper, we shall consider helium 4 un- 
less it is explicitly specified that it is helium 3. In the next 
sections, we show that shock waves can occur in this sys- 
tem, and we compute the radius (denoted as shock length 
in the foUowings) at which a shock wave occurs for various 
oscillation amplitudes Axq of the transducer surface. 



2.2 The method of characteristics 

In the study of compressible fluids, the method of charac- 
teristics is a standard one It has already been used 
for helium in a planar geometry (see for example |6)). 

Here we are interested in the shock length in spherical 
geometry, for the equation of state (|l|). In this paper we 
define the shock length as the distance from the center 
where the shock forms. We shall use the method of char- 
acteristics to predict a lower bound for the shock length. 



2.2.1 Rieman invariants and characteristics 

Let us first recall the principle of the method. Solving the 
Euler equations means knowing the density p and the ve- 
locity u everywhere within a certain domain of the (r,t) 
plane. In our case, the "characteristics" are two families 




Fig. 1. The two characteristic families defined in the {r,t) 
plane. Here each curve is labeled by the r-coordinate of its 
intersection with the t = axis. The point indicated by a 
small circle can be referred to either by its coordinates (r,t), 
or by the labels of the characteristics which intersect at this 
point (i,j). 



of curves jc,-^"* | and Ic,,- H parameterized by i, 

where the parameter i could for example be defined as 
the r-coordinate at which the characteristic cuts the axis 
t — (see Fig. || for an illustration). Thus, C-^^ (resp. 

c| •*) refers to a single curve in the (r,t) plane, which be- 
longs to the (+) family (resp. (— )), and has in our case 
a positive slope (resp. negative). Each family completely 
covers the (r,t) plane as i is varied. The two families of 
curves intersect each other. Then, instead of locating a 
point in the plane by its coordinates (r,t), it is equivalent 

to indicate which particular characteristics C^-^^ and 
intersect each other at this point. The point (r,t) will then 
equivalently be referred to as point («, j). 

These families are chosen so that, for each family |c|'''' | 

with fc = ±, there exists a quantity Ik called "Rieman 
invariant" , which is a function of p and u, and obeys a 

(k) 

simple evolution equation along any characteristic Q of 
the family. If the value of Ik is known at one point of a 
characteristic Cp'' (for example at the initial time), then 
it is easy to compute it on the whole curve. 

As each point of the (r,t) plane is the intersection of 



two characteristics C^^-* and Cj '', we know the values of 
I+{p, u) and /-(p, u) at this point. Then the density p and 
the velocity u are entirely determined. 

Now the precise form of the Rieman invariants and 
corresponding characteristics have to be derived from the 
Euler equations (^) and the equation of state for helium 
(^. The detailed calculations are given in Appendix A. 
Here we just give the results concerning the shape of the 
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characteristics and the Rieman invariants. By definition, 
at any point, the slope of a characteristic C^^"* with equa- 
tion r = x{t) is dx/dt = u + Cs where u and Cg are taken 

in (x(t),t). For a characteristic c\ \ it is dx/dt = u — c^. 
The sound speed Cg is given from the equation of state 

by 

cs = a(p-i). (7) 

Then the derivative along the characteristic reads 



— = 5t + (-U ± Cs)dr. 

The Rieman invariants are found to be 
I± = Co (p - In p) ± u 



(8) 



(9) 



The equations verified by the "invariants" Ik with k = 
± read 

d , , C„(p — l)u 

-[i,]^e^^^^o (10) 

where the derivative ^ is taken along a characteristic c| . 
Again, 9 is respectively 0, 1, or 2 in planar, cylindrical, or 
spherical geometry. 

In the case of planar geometry [9 — 0), Ik is a true in- 
variant, since it is constant along the corresponding char- 
acteristic, hence its name. In spherical or cylindrical ge- 
ometries, it is not constant, due to the source term in 
Eq.Q, though it is still called an "invariant" . 

All this is valid at least as long as the characteristics 
belonging to the same family do not intersect each other. 
When they do, the corresponding Ik becomes multival- 
ued. This is an indication that a shock has formed, and 
beyond the corresponding time, the description used in 
this section breaks down. 



2.2.2 Lower bound for the shock length - Analytic 
calculation 

We shall now show that such an intersection does occur in 
the system, and calculate the corresponding shock length. 
We consider a spherical domain bounded by a spherical 
piston. At t = 0, the fluid is at rest with a density pst (label 
"st" standing for static) and the piston surface is a sphere 
of radius Lq. As long as the fluid is at rest, all the char- 
acteristics are straight lines (see Fig. |^). The respective 

RTnd|cf"^| are 

-|-Cst and — Cst, Cst being the sound velocity for the initial 
density pst- The piston starts to move at t = with a 
velocity Vp{t) = —Avosm{ojt). 



slopes of the characteristics j C]'^' )■ and 



We denote Cg ^ the characteristic which originates from 
= Lq when t=0, with slope —Cst- The domain to the 



left of Cq ^ is unperturbed unless some characteristic C- 

crosses Cq \ i.e. unless there is a shock. 

The aim of this calculation is to find an upper bound 
for the time necessary to form a first shock in the system. 
As the piston moves, it emits some characteristics which 
will cut Cg after a while, leading to a shock. We only 



(-) 



0.20 



0.15 



0.10 



0.05 



0.00 




Fig. 2. Characteristics obtained for an experimental oscilla- 
tion amplitude equal to Axq = 10 /i m, a cell radius Lq/X = 10, 
and a time step St — 0.01 T. The location of the piston is rep- 
resented by the solid thick line on the right. On the upper left, 
two characteristics are crossing each other and the program 
stops. We chose a much greater oscillation than in the exper- 
iment, in order to have a rapid shock formation, and thus a 
readable picture. 



study how the characteristics emitted at early times, and 
almost parallel to Cg \ will eventually cross it. Of course, 

some other characteristics emitted later could cross Cg"'' 
earlier, or some shocks could occur somewhere else at ear- 
lier times. That is why our calculation only gives an upper 
bound for the shock time. 

The details of the calculation are given in Appendix B. 
We find that an upper bound for the time of shock forma- 
tion is 



hock 



< 



Cst 



1 — exp 



Pst 



1 



2 Lo cj Avo Pst 



< 



Lo 



Cst 

As the corresponding shock length rghock is measured from 
the center of the sphere, a lower bound for rshock is 



Lo > ^shock > Lq exp 



Pst - 1 



2 Lo uj Avo Pst 



> 0. (12) 



From this study we conclude that in spherical geom- 
etry, there is always formation of a shock, whatever the 
velocity of the piston is, as long as it has a nonzero accel- 
eration towards the center of the sphere. However, when 
the oscillation amplitude Axq = Avo/ij-i goes to zero, 
the shock forms very near the focal point. As shocks are 
formed by the intersection of tangential characteristics, 
their initial amplitude is zero. If they are formed very near 
the center, their amplitude does not have time to grow 
much. This is true especially because of the existence of 
a cut-off: the notion of shock becomes meaningless when 
the width of the shock becomes of the same size as the 
shock length itself. So, as the oscillation amplitude tends 
to zero, the jump in density and velocity at the shock 
also vanishes, and their is no contradiction with the fact 
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Fig. 3. Schematic representation of the network formed by the 
two families of characteristics, where i and j are integers. 

that the solution is expected to approach the hnear solu- 
tion 

The above calculation is analytic and gives only a lower 
bound for the shock length because we have assumed that 
the shock occurs on the first characteristic Cq ^ 

2.2.3 Numerical calculation of the shock length 
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Fig. 4. Shock distance in helium 4 versus the amplitude of the 
oscillation on the ceramic, for a cell radius Lo — 33. 6A. We 
compare the analytic prediction (solid line) and our numer- 
ical calculations (symbols), both based on the characteristics 
method. The inset shows an expanded view for small ampli- 
tudes. The arrow indicates the amplitude for which cavitation 
is observed in the experiment, starting from a zero static pres- 
sure. 



As we shall see now, a numerical simulation shows that 
our calculation of the lower bound is in fact a good esti- 
mate of the shock length itself. Let us now solve Euler's 
equations numerically by computing a network of char- 
acteristics |c,-~H and ■!cf'''H . The values of the 
Rieman invariants, and thus of the density p and velocity 
u, will be defined on the intersections of the left cj "* and 

right cf^'' characteristics. As the fiuid is at rest on the left 
of Cq \ we restrict our calculation to the {r,t) domain 
located between Cq and the piston trajectory. 

Several parameterization choices are possible. For the 
family, we have chosen to take i — ti/6t for each 

characteristic cl"*"^ , where ti is the time at which C^^' and 

Cq intersect, and St is an arbitrary fixed time step. We 
discretize the system by restricting i to integer values. 
This means that the subset of characteristics that will be 
considered for the simulation are initialized at regular time 
intervals iSt when they cross the first characteristic Cq 
emitted by the piston (see Fig. ||). 

When C^^-* meets the piston, a new characteristic C^- ■* 
(with the same label i) is emitted from the piston at the 
same time. This defines the parameterization of the (— ) 
family. The intersections of these two families of charac- 
teristics form a network of points whose locations will be 
determined in the following. The intersection of the char- 
acteristics c|^' and cj""* is referred to as point 

During the n*'' step, we compute all the points (?,j) 
such that i -\- j = n with i and j integers. Let us describe 
now how a point {i, j) can be computed from the points of 
the previous step (see Fig. ||). All the information comes 



from the two sites — 1) and (i — We must ex- 

trapolate each characteristic C,^^'' and cj up to the next 
intersection (i, j). We compute the local slopes u± Cs of 
the characteristics C^^^ and cj •* in sites — 1) and 
(i — 1, j) respectively. Then («, j) is the intersection of the 
two straight lines which respectively go through (i, j — 1) 
and (i — 1, j) with these slopes. 

The values of the Rieman invariants /+ and /_ at (z, ,7) 
are found by numerical integration of the equations (pUlQ) 
along the two involved characteristics. From these values, 
both p and j = pu can be obtained. 

The program stops whenever two characteristics of the 
same family cross each other, as shown on Fig. |^. Then a 
shock occurs and the calculation based on characteristics 
breaks down (there would be multivalued points). 

Let us now compare analytic and numerical results. 
Both calculations were done by taking Xq equal to the ex- 
perimental cell radius, i.e. 8 mm = 33.6 A. In Fig. ^, we 
plot the shock distance rghock measured from the center of 
the sphere as a function of the amplitude of the oscillation 
of the piston. We compare it with the analytic result. The 
agreement is excellent. This shows that our analytic cal- 
culation not only gives a lower bound but in fact a good 
estimate for the shock distance itself. 

We have also performed the analytic calculation in the 
case of helium 3, neglecting the role of viscosity. It is ex- 
pected that shocks form at smaller wave amplitude in he- 
lium 3, because, in this lighter liquid, the spinodal pressure 
is less negative than in helium 4. 
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Fig. 5. Staggered lattice used for the numerical simulations. 
The radius of the simulation domain is L — (if + ^) Sr. The 
symmetry with respect to the center r = imposes jo ~ — J-i, 
so that the velocity would vanish at the center : j(r=0) — 0. 



2.3 Numerical simulations of Euler's equation 

2.3.1 Numerical method 

As we shall see now, a simple finite difference numerical 
scheme is sufficient to simulate our system with moderate 
amplitude, at least as long as the shock length is smaller 
than the spatial discretisation step. Our aim is to calcu- 
late numerically the pressure and density oscillation at 
the center. We chose to have two lattices, one for mass 
and the other for momentum. They are staggered (Fig. ||) 
and allow us to enforce exactly the conservation of mass: 



t+5t/2 
Pk 



t-St/2 
Pk 



5t 



Sr 



j_ 

12 



We have taken the momentum equation in the form: 



(13) 



■t+st 

Jk 



fk 



^0 



5t 
Sr 



. t+St/2 



, t+&t/2 
^Pk 



' t+St/2 
Pk 



t+5t/2 
Pk+1 



- 1 



Pk 



t+St/2 



t+5t/2 
Pk 



Sr 



St 



St 



Sr 



t+5t/2 
Pk 



t+5t/2 
■ Pk+1 



(14) 



As we use a staggered lattice, we only have to specify the 
boundary conditions for the momentum. It is vanishing 
at the center of the sphere (r=0), and thus the symme- 
try with respect to the center imposes jo — —j-i (see 
Fig. H). On the other hand, the motion of the piston is 
implemented by 

.jK{t) = —pK^'^XQ sin([x)t) [1 — exp(— i/1.5)] . 



2.3.2 Focal pressure 

In Fig. ^, the focal pressure is represented as a function 
of time. It is calculated from the average density in the 
central cell. The results of Fig. |6| were obtained for a cell 
length equal to the experimental one, i.e. Lq = 33.6 A, 
an oscillation amplitude Axo = 0.7 nm, and a zero static 
pressure. Then the reduced sound velocity is equal to 1, 




t/T 

Fig. 6. Focal pressure for a weak oscillation Axq — 0.7 nm. 
The simulation was done with 100 mesh points per wavelength, 
starting with a static density pst = po, and thus a vanishing 
static pressure. 




Fig. 7. Same figure as Fig. |^ with a zoom on the steady state 
region. We are nearly m the linear regime. 



and the wave needs 33.6 time units to reach the center of 
the cell. We are interested in the steady regime, which is 
established around t/T > 45. For t/T > 33.6, the wave re- 
fiected on the focal point propagates towards the ceramic, 
and then back to the focal point. It reaches the latter at 
t/T = 100.8 = 3 X 33.6. Then we stop our measurements, 
to be consistent with the experiment, which uses short 
bursts. 

As long as the oscillation amplitude of the ceramic is 
not too large, non-linear effects are negligible far from the 
focal point. If we observe the focal pressure during the 
steady regime, we find it also nearly sinusoidal, and the 
positive swings are only slightly larger than the negative 
ones (Fig. 0). 

We have checked that the amplitude of the density is 
well represented by the function sin(A;r)/(fcr), as predicted 
from the linear theory jlj] (see Fig. ||). 
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Fig. 8. Density profile for a weak oscillation Axq — 0.7 nm 
(same simulation as in Fig. |^. We are still in the linear 
regime. The thin solid lines indicate a fit of the amplitude by 
1/r and the dashed line a fit by sm{kr) / (kr) . The latter is 
almost identical to the numerical result, except near the focal 
point where non-hnear effects become visible. The simulation 
was done with 100 mesh points per wavelength, starting with a 
static density pst = po - and thus a vanishing static pressure. 




From our calculation of section p.2| , we expect a shock 
to occur near the focal point. However, for the case of 
Fig. ||, it happens on a region around r — that cannot be 
seen because it is much smaller than the mesh size. Indeed, 
the predicted reduced shock length would be Tc/A = 4.2 x 
10"^^, to be compared with the mesh size 6r/X = 0.01. 

When Axq is increased to 5 nm, non-linear effects be- 
come more important. One sees the formation of fronts 
(Fig. H). The reduced shock length (from the center) is 
equal to 0.003 and is thus still smaller than the spatial 
step 6r/X = 10~^. At the center, the positive swings of 
the pressure are now much larger than the negative ones 
(Fig. |l^). There are several sources for non linearities: 
the equation of state; the inertial term in Euler equations; 
both enhanced by the spherical geometry. The relative im- 
portance of these factors will be discussed now. We have 
performed various simulations where we suppressed one 
term or the other. Results are summarized in Table |^. 

If we take a constant sound speed Cg = Cgt, i.e. a linear 
equation of state, then non-linear effects are strongly re- 
duced and the maximal and minimal pressure excursions 
are almost symmetrical. 

If we rather suppress the inertial term u ■ Vw from the 
Euler equation for the velocity, non-linear effects are also 
reduced, but to a lesser extent. 

Thus all non-linear terms reinforce each other, but the 
dominant non-linear effect comes from the equation of 
state (the geometry also plays a crucial role of course, 
which cannot be separated from the other effects). 

Thi s is co nfirmed if we do the analytic calculation of 
section 2.2.2 again, now with a constant sound speed. We 




Fig. 9. Density profile for an oscillation Axq — 5 nm, at dif- 
ferent times, corresponding to the maximal and minimal focal 
pressure. Non linearities are becoming important. The simula- 
tion was done with 100 mesh points per wavelength, starting 
with a static density pst = po- 




t/T 

Fig. 10. Focal pressure for an oscillation Axq — 5 nm. Non 

linearities are becoming important. The simulation was done 
with 100 mesh points per wavelength, starting with a static den- 
sity Pst = po ■ 



Lo 1 Ao 


Lo (mm) 




case 1 


case 2 


case 3 


10 


2.38 


Pmax 


13.53 


6.18 


11.92 






J mm 


-3.49 


-4.83 


-3.58 


20 


4.76 


P 

max 


17.17 


6.34 


14.27 






J mm 


-3.31 


-4.73 


-3.41 



Table 1. Maxima and minima of the focal pressure (bars) 
computed in three different cases. Case 1: full simulation of 
the equations (^). Case 2: the sound speed is kept constant 
Cs ~ Cat. Case 3: the inertial term u ■ Vit is suppressed. All 
calculations are done for the same experimental oscillation am- 
plitude Axo = 5.9 nm, and for two different cell radii. 
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Fig. 11. Shock distance in helium 4 versus the amplitude of 
the oscillation of the ceramic. We compare the analytic predic- 
tion for the full equation of state (solid line) and the case of a 
constant sound speed (dashed line). 



find then that 



J^shock ^ Lq exp 



LoujAvq 



(15) 



On F ig. 11, we compare with the resuh ( [l2[ ) of section 
2.2.2 obtained with the full equation of state. It turns out 
that, for reahstic oscillation amplitudes, the shock appears 
almost at the center if the nonlinearities of the equation 
of state are not taken into account. 



2.4 Shock formation at large amplitudes. The WENO 
scheme 

Eventually, when the amplitude Axq is further increased, 
the formation of shocks can be observed in our simula- 
tions. The numerical scheme described in the previous 
section becomes unstable when rshock — Sr, i.e. for a van- 
ishing static pressure and a reduced Sr = 10^^, when 
Axq — 6 nm. Then a new numerical scheme has to be 
used. 

For performing fine analysis of the Euler flow dynam- 
ics, the numerical scheme must recover low dissipative pro- 
perty. Non-dissipative high-order accurate schemes (like 
spectral or Fade schemes) have been identified as suitable 
tools as far as regular numerical solutions are searched. 
Nevertheless, when dealing with compressible flows involv- 
ing discontinuities, non-dissipative high-order schemes in- 
troduce spurious oscillations in the vicinity of the dis- 
continuity and one must use a numerical scheme which 
can both represent the smooth regions of the solution 
with the minimum of numerical dissipation, and capture 
the discontinuities by using an ad hoc scheme with a ro- 
bust discontinuity-capturing features. Therefore, as shock 
waves may occur in the computational domain for the 
present calculations, the numerical method we use is based 



on a high-order intrinsically-dissipative scheme originally 
designed to capture discontinuities. The method, called 
WENO (Weighted Essentially Non-OsciUatory) @, is pre- 
sented in the next section, and some more technical details 
about WENO and ENO schemes can be found in the lit- 
erature @,^|l|,0,|g 

As in many schemes devoted to computations involving 
shocks, the WENO scheme uses the Riemann invariants 
as variables. Computing the evolution of these variables 
requires to know their values not only at integer space 
coordinates, but also at some intermediate locations. An 
extrapolation from integer positions is thus necessary and 
this is where the fundamental idea of WENO schemes 
comes in. The time evolution scheme (a third order Runge- 
Kutta) and the change of variables towards the character- 
istic plane are more classical, but we shall still recall them 
for self-consistency. 



2.4.1 The WENO method in more details 

For simplicity, the governing equations (|3|) are recast in 
the following abridged form: 



dQ 

at 



C (r, Q) 



(16) 



where Q — {p , p u) is the vector of the conservative vari- 
ables and 

or 

stands for a spatial operator, applied on Q, based on both 
the Euler flux vector F (Q) — (^p u , p + and the 

source term S (r, Q) = — (p u, p u^)*. 

In view of the discretization of the Euler equations 
(|l6|), we will denote by St and Sr the time step and cell 
width respectively. Q" will denote the numerical vector 
solution at a time t = to + n ■ St and at a position r = 
i ■ Sr. For simplicity, the integration of these equations has 
been performed by means of a decoupled time and space 
algorithm. 

a) Time integration 

The time integration is then performed by means of 
a third-order accurate Runge-Kutta method, proposed by 
Shu and Osher iQ, chosen because this high order ac- 
curate scheme does not increase the Total Variation of 
the right-hand-side of the equations {£ (r, Q)). When deal- 
ing with discontinuities, this Total Variation Diminishing 
(TVD) property is important since it ensures that no lo- 
cal extremum can be created during the time integration, 
meaning that non oscillation may occur in the shock wave 
vicinity due to the time scheme. At each point of the com- 
putational grid, the time integration is then obtained via 
a multi-step algorithm as follows: 



Q 



(0) 
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(0) 



(0) 



St C 
1 

+ 4 
2 



i-2 



i-1 



i+2 



i+3 



(5< /:(q'^ 



(17) 



-dt c 



The explicit Runge-Kutta scheme exhibits a stability 
condition based on the CFL number: CFL = fj, 6t/6r, 
where is the maximal value of the eigenvalues ^ J , to be 
defined below. All the calculations presented herein have 
been obtained by considering CFL = 0.5, which corre- 
sponds to a nearly optimal value for the considered Runge- 
Kutta scheme. Let us mention that this value ensures a 
good representation of all the time scales of the flow con- 
sidered. 

b) Spatial integration 

The spatial discretization of the right-hand-side term 
C (r, Q) of equations (|l6|) is obtained by means of a high- 
order finite difference scheme: 



1 

Sr 



F 



-1/2 



F 



-1/2 



Sir.,Q7) (18) 



where is the numerical flux evaluated at the cell 

interface (?'i+i/2)- To reconstruct the numerical flux at the 
cell interfaces, a scheme with a discontinuity-capturing 
feature must be employed to prevent oscillation in the 
vicinity of the shock wave. Following a previous study |0 
on the capability of some recent high-order shock captur- 
ing schemes to recover basic fluid mechanic phenomena, 
the numerical flux has been evaluated by means of Essen- 
tially Non-Oscillatory (ENO) family scheme [0|l|,|l3. 
The numerical flux is approximated by means of poly- 
nomial reconstruction over several grid points (the set of 
these points is named "stencil") around the cell interface 
considered. We shall now describe this reconstruction of 
the fluxes in details. 

• Change of variable: For simplicity and accuracy 
purposes, the discretization of the Euler flux is based on 
a polynomial reconstruction applied on the l ocal c harac- 
teristic variables (Riemann invariants, see § 2.2.1) since 
the equations recover the scalar form. Then, the propaga- 
tion directions can easily be followed in the characteris- 
tic plane. In order to perform this change of variables, 
the method is the same as in Appendix A. First, one 
linearizes the Euler equations and finds the eigenvectors 
of the Euler flux jacobian evaluated at the cell interface 
{{dF/dQ)^j^^l2)- l^ote that, as we are now using the Eu- 
ler equations in conservative form (^) instead of (||), the 
jacobian differs from the one given in appendix A (matrix 
A), and reads now 







1 



dF/dQ = 



-J_ 

n2 



(19) 



In order to compute {dF/dQ)-_^^^,2- - and then the eigen- 



values (/i; 



4 + 1/2^ 



and the left i'^i_^_l/2) and right (r*^^^ 



i+1/2 



o 



o 
o 



o 
o 



o 

© o 
©00 



} 



|i>0 



} 



|i<0 



Fig. 12. Sketch of all the stencil candidates to recover a third- 
order reconstruction (p — 3) at the cell interface ri+i/2- The 
eigenvalue /i stands for the extrapolated value fJ-^^i 



+ 1/2- 



eigenvectors -, the conservative variables Qf must be eval- 
uated at the cell interface rj_|_i/2 (note that a bold 
refers to the eigenvector, while ri+1/2 stands for the scalar 
spatial coordinate). As these variables do not vary linearly 
in the cell if a shock is present, one cannot use a simple 
arithmetic or geometric average, but rather a Roe average, 
whose description can be found in Ref. |]l^. This ensures 
the consistency of the scheme, i.e. that {dF/dQ)^^-j^/2 con- 
verges towards {dF/dQ)^ when Sr tends to zero. 

The numerical Euler fiux is then projected onto the left 



eigenvector matrix (lj^|/2' "^^^ scalar ENO recon- 

struction procedure is applied to the projected fluxes [^6|. 
In the physical domain, the numerical Euler flux are then 
obtained by a projection onto the right eigenvectors 



(-) 



-1/2 



^2+^1/2) ^^"^ read: 



F, 



+1/2 



k=± 



ENO" 



1+1/2 • '■i+1/2 



(20) 



1/2) 



where f^^'~' i+1/2 stands for the scalar ENO reconstruc- 
tion, which will be defined below. 

• Extrapolation of variables at non-integer lo- 
cations - the core of (W)ENO schemes: 

First, we shall present the strategy used in ENO (Es- 
sentially Non-Oscillatory) schemes |16|. For a given non 
integer location, there are several ways to perform the ex- 
trapolation, depending on from how many integer points 
it will be performed, and how these points will be located 
with respect to the non- integer one. The set of these points 
is called a stencil. Fig. |l2| illustrates the different possible 
choices for a stencil of length p = 3. Simple finite differ- 
ence schemes use a stencil defined once for all. Here, only 
the length p of the stencil is fixed for a given simulation. 
All the p + 1 possible locations are considered as candi- 
dates, provided that there is at least one point adjacent 
to the non-integer point (see again Fig. [l^ ). 

If we denote by p the order of the reconstruction, the 
ENO procedure |l^ chooses the most regular stencil among 
the p + 1 stencil candidates. As an exemple for p = 3, we 
can see all the stencil candidates on Fig. |l^. A first se- 
lection among the p + 1 stencil candidates is performed 
according to the sign of the two eigenvalues (Ai*'Yi/2)- '^'^^ 
keeps the p most left stencils for the positive or null eigen- 
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p 


m 


j = 


3 = 1 


j = 2 


j = 3, 


j = 4 


1 


-1 


1 















1 










2 


-1 


-1/2 


3/2 













1/2 


1/2 










+ 1 


3/2 


-1/2 








3 


-2 


1/3 


-7/6 


11/6 








-1 


-1/6 


5/6 


1/3 











1/3 


5/6 


-1/6 








+ 1 


11/6 


-7/6 


1/3 






4 


-3 


-1/4 


13/12 


-23/12 


25/12 






-2 


1/12 


-5/12 


13/12 


1/4 






-1 


-1/12 


7/12 


7/12 


-1/12 









1/4 


13/12 


-5/12 


1/12 






+ 1 


25/12 


-23/12 


13/12 


-1/4 




5 


-4 


1/5 


-21/20 


137/60 


-163/60 


137/60 




-3 


-1/20 


17/60 


-43/60 


77/60 


1/5 




-2 


1/30 


-13/60 


47/60 


9/20 


-1/20 




-1 


-1/20 


9/20 


47/60 


-13/60 


1/30 







1/5 


77/60 


-43/60 


17/60 


-1/20 




+ 1 


137/60 


-163/60 


137/60 


-21/20 


1/5 



Table 2. The constant coefficients of the ENO reconstruc- 
tion up to 5*''-order. 



values and the p most right otherwise (Fig. Indeed, 
the sign of the eigenvalue gives the propagation direction 
of the associated characteristics, and thus of the relevant 
information. 

The regularity of the function on each of the p re- 
maining stencils is measured by the undivided difference 
Table ijl^l evaluated on each stencil and the most regular 
stencil is chosen among all the p stencil candidates. Of 
course, this stencil may be different at each time step, for 
each location, and for each eigenvalue (m*'Yi/2)'==±- 

The scalar ENO reconstruction is then applied on this 
specific stencil by means of the following polynomial de- 
velopment: 



cENO 



i+l/2 = X! ■ ^ (Qi+m+j) 



(21) 



The integer m refers to the index of the most left point 
of the chosen stencil. The sum goes over all the points 
of this stencil. One recognizes the projection of the fluxes 
onto the left eigenvectors, that allows to go from physical 
variables to characteristic variables. The constant coeffi- 
cients of the polynomial (Cj'*^) are calculated in order to 
recover a scheme of order p in regular regions. The values 
of Cj can be found in Table || up to p = 5. 

Using this generic ENO scheme (^, 21), the Eulcr flux 
derivative is estimated with ap*''-order of accuracy at best 
(in regular regions). However, when the stencil used at 
the cell interface ?'j-(-i/2 is different from the one at rj_]^/2 
(which is the case in strong gradients or shock regions), 
the order of accuracy decreases. 

One of the drawbacks in the generic ENO scheme is 
the necessity to check and choose between p stencil can- 



p 


n = 


n = 1 


n = 2 


n — 3 


71 = 4 


2 


1/3 


2/3 








3 


1/10 


6/10 


3/10 






4 


1/35 


12/35 


18/35 


4/35 




5 


1/126 


10/63 


10/21 


20/63 


5/126 



Table 3. The constant coefficients of the WENO recon- 
struction up to 5*''-order, for a positive eigenvalue. 



didate, which is quite CPU consuming. To overcome this 
disadvantage, we preferred using a Weighted ENO scheme 
(WENO) since it improves the order of accuracy of the 
generic ENO scheme by using a weighted combination of 
the p possible stencils. The weights depend on the 
degree of regularity of the solution. In regular regions, 
they can be computed to achieve {2p— l)*'*-order of accu- 
racy whereas in regions with discontinuities they are set 
to zero, leading to a standard ENO scheme. The WENO 
flux is estimated by: 



cWENO'^ 



i+l/2 



p-1 
n=0 



rENO 



k,n—p+l 
+ 1/2 



(22) 



where f^^'^\'-^i/2^^ is the generic ENO reconstruction 



given by equation 
as follows: 



11) and ojn are the weights defined 



n=0 



p-1 
1=0 



with Pn = 



CP 



(e + ISnY 



(23) 
(24) 



e is a small positive number to avoid denominator to be 
zero (hereafter we set e = 10~^) and ISn is a measure of 
the flux function regularity for the rt"* ENO stencil can- 
didates. The evaluation of the smoothness measurement 
{ISn) is based on the undivided-differences and the ISn 
formulation can be found in . It is such that a more reg- 
ular curve gives a smallest ISn , and thus a largest weight 
Pn- The CP coefficients are reported in Table |^ up to the 
order p — 5, for positive eigenvalues at the cell interface 
(^^i+l/2 *^)- ^'^^ negative eigenvalue case, the WENO 
coefficients can be obtained by symmetry with respect to 
the considered cell interface (?Y+i/2)- 

Using this WENO scheme (|2^EJ) , the Euler flux deriva- 
tive is estimated with a (2p — 1) -order of accuracy at 
best (in regular regions). Moreover, let us underline that, 
if the solution is regular enough, the WENO procedure 
recovers a high-order centered scheme, which is of course 
non-dissipative . 

If the length of the stencil is p = 3 - which is our 
case, the order of the scheme in regular regions is then 5. 
However, as the order of the scheme drops to 1 in shock 
regions - as for any other scheme - it is not interesting 
to increase too much the order in regular regions. As the 
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scheme is less dissipative with p = 4 in the non regular 
regions, this could even make the results worse (indeed, in 
our case, p — A gave less satisfying results), since spurious 
oscillations may occur. 

c) Boundary conditions 

To solve the system of equations we need boundary 
conditions at the center of the sphere (r = 0) and at the 
ceramic surface (r — L{t)). At the center of the sphere, 
since the velocity of the fluid is an antisymmetric quan- 
tity, the fluid is at rest (wl^^g ~ Moreover, as the pres- 
sure and the density are symmetric quantities, their radial 

d 

derivatives are set to zero ( — (P, p) 



2.0 



1.9 



—1 1 1 r- 



dr 



0). The singu- 



r=0 



lar behavior of the source terms involves a specific treat- 
ment at the sphere center. While the velocity tends to 

pVb 

zero when the radius tends to zero, the ratio — tends 

r 
pu 

to a limit which has to be determined though tends 

r 

to zero. As — is a symmetric quantity, its limit (reached 
r 

at r = 0) is calculated by writing that its radial deriva- 
tive vanishes at the sphere center by using a second order 
upwind difference: 



1.4 




2 

r/l 



Fig. 13. Density profile for an oscillation Axq = 7 nm, at dif- 
ferent times, corresponding to the maximal (solid line), mini- 
mal (dot dashed line) focal pressure, and to a sharp front ar- 
riving to the center (dashed line) just before the pressure maxi- 
mum. Non linearities are becoming very important. The simu- 
lation was done with 100 mesh points per wavelength, starting 
with a static density pst = po- 



d^pu. 
or \ r / 



_ 1 
= 0. 



(v)«-(?)„.-( 



'pu\ 
r J i=2 



(25) 



At the ceramic surface, one has to impose a condition 
corresponding to the motion of the ceramics. In our case, 
as the displacement of the sphere does not reach a sonic 
velocity, the two eigenvalues and p~ are of opposite 
sign. This means that the two informations necessary to 
determine the two unknowns p and u come from opposite 
directions. In particular, at the boundary, one comes from 
the interior of the domain and the other from the outside. 
Thus, it is possible to prescribe one of the variable, and 
the other one will be determined by an upwind scheme, 
i.e. a scheme asymmetric towards the inner domain. In the 
present problem, it is much more natural to prescribe the 
velocity at the ceramic surface: 

u {L{t),t) = -u Axq sin(cjt) [1 - exp(-t/1.5)] . (26) 

Axq is the amplitude of the displacement of the ceramic 
and Lo is the ceramic pulsation. These two parameters are, 
of course, prescribed. The exponential ter m re presents the 
response time of the transducer (see Sec. 2.1). 




t/T 

Fig. 14. Focal pressure for an oscillation Axq = 7 nm. Non 

linearities are becoming very important. The simulation was 
done with 100 mesh points per wavelength, starting with a static 

' Pst = Po • 



2.4.2 Numerical results 



Of course, all the results presented in section 2.3.2 can be 
reproduced with our WENO scheme. Besides, one can in- 
crease further the oscillation amplitude. For Axq = 7 nm, 
the reduced shock length (from the center) is equal to 
0.043, which is now larger than the spatial step 5r/X — 
10^^. One clearly sees the formation of sharp fronts (Fig. ^ 
At the center, the positive swings of the pressure are now 



very sharp compared to the negative ones (Fig. Q). The 
shape of the negative swings is also clearly asymmetric in 
time. 

It would be tempting, in order to reduce the compu- 
tational effort, to assume that non-linearities play a role 
only in the last wavelengths. Then, one could simulate a 
reduced box with radius Lrcd, and take as an input con- 
dition : 

Axrcd = Axcxp-;^ (27) 



^rcd 



, where Axcxp is the experimental oscillation amplitude of 
the transducer, and Lcxp its radius. Actually, this is fine as 
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Fig. 15. Minimal and maximal focal pressure, as a function 
of the cell radius. The product Lrcd x Ax^ad is constant for 
each curve. Solid lines (resp. dashed lines) correspond to an 
experimental oscillation amplitude zixcxp = 0.7 nm (resp. 7 
nm). One of the curves (Pmax for 7 nm) is not visible, as it 
varies from 28000 to 106 bars! In all cases, the simulation was 
done with about 350 points per wavelength. 



long as non-linearities do not play a role at all, i.e. as long 
as the focal pressure is sinusoidal. In all other cases, as 
this is shown on Fig. non-linearities are built through 
the whole propagation process, and one cannot neglect 
them even far from the center, without modifying the fo- 
cal pressure. This could also be seen in the shock length 

which is directly proportional to 



analytic expression 
the cell radius. In Fig. 15|, we have performed simulations 
with various cell radii, corresponding to the same exper- 
imental oscillation amplitude on the ceramic. As the cell 
radius increases, both the positive and negative pressure 
swings decrease (in absolute values) . Non-linearities make 
it more difficult to reach extreme values. Thus in all the 
simulations presented in this paper, we have simulated the 
whole experimental cell, with radius 8 mm. 

We have seen on Fig. |lj that positive pressure peaks 
can reach tremendously high values. Of course one may 
wonder whether this is physical. The obvious answer is no 
- but it is worth to discuss this point in some details. 
The solution of the problem as we defined it in Section 



2.1 



becomes singular when the shock reaches the center of 
the sphere, and positive peaks actually diverge in the sim- 
ulations as the spatial discretisation step 5r is decreased 
(see Fig. |l^). For each fixed 5r, one still finds a finite focal 
pressure, as it is defined as an average over the central cell 
- with radius 5r/2. We check on Fig. |l^ that the minimal 
pressure does converge for a decreasing Sr. 

A first remark is that this singularity involves only a 
very small region around r = (ultimately, it is singular 
only in r = 0). In the experiment, one cannot measure the 
pressure exactly in r = 0, but rather over the whole region 
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Fig. 16. Maximum of the focal pressure (solid line) as a func- 
tion of the number of points per wavelength, for an oscillation 
amplitude Axq — 7.0 nm. For comparison, the dashed line gives 
the pressure after a spatial average weighted by a Gaussian with 
waist 7/im. 
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Fig. 17. Minimum of the focal pressure (solid line) as a func- 
tion of the number of points per wavelength, for an oscillation 
amplitude Axo — 7.0 nm. For comparison, the dashed line gives 
the pressure after a spatial average weighted by a Gaussian with 
waist 7/im. 



reached by the laser beam (see section |[). The intensity of 
the laser beam through its cross section is Gaussian, with 
a waist (half width) equal to 7 microns. In order to take 
this averaging effect into account in the simulation, we 
have also computed a spatial average value of the pressure 
weighted by a Gaussian of waist Tfim (Fig. |8|). Positive 
pressure peaks are lowered, while the negative swing is not 
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t/T 

Fig. 18. Averaged focal pressure for an oscillation Axo — 7 
nm. The spatial average is weighted by a Gaussian with waist 
7fim. This is to be compared with Fig. \lA 



very sensitive to the averaging process, as gradients are 
much weaker than during the positive swing. But above 
all, one now obtains extrema which are converging as Sr — > 
(Fig. p. 

But it is not enough to rule out the singularity by an 
averaging process. Actually, the singularity is not expected 
to hold as such in a more realistic description. First, the 
third Euler equation for energy should also be taken into 
account in this regime, as well as regularization mecha- 
nisms (dispersion, dissipation). Secondly, in the experi- 
ment, we expect diffraction, and the fact that actually the 
flux is not zero at the focal point, to break the geometrical 
symmetry. Making quantitative estimates of these effects 
is not simple, and is postponed to future work. Finally, 
even with averaging, positive pressure peaks in Fig. 08^ 
are still so high that they are far above the solidification 
pressure (25.3 bar at T=0). As we shall see in the last 
section, our recent experiments show that indeed acoustic 
waves can trigger crystallization, not only cavitation Q. 
In our simulations, we are not able to treat the possible 
crystallization. 

As a conclusion, the value of the positive pressure 
swings is not expected to be reliable for high amplitudes 
of the transducer surface. But we checked, using different 
types of numerical regularization of the shocks (an exam- 
ple is given in the next section), that even if it affects the 
maximal pressure, it has no effect on the negative swings, 
and thus it is still possible to use our simulations to draw 
conclusions for negative pressures. 

In Figs. |l^ and |2^, we show the density profile and fo- 
cal pressure obtained for Axq — 30 nm, a value which can 
be reached in cavitation experiments. Fig. 19 illustrates 
how shocks are formed when the wave arrives near the fo- 
cal point. The amplitude of the shock increases when the 
shock itself arrives at the focal point, leading to tremen- 
dous pressure maxima in the simulations (here about 12000 
bars). This is of course unphysical, as we just discussed it, 
but the important point is that the minimum (negative) 
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Fig. 19. Density profile for an oscillation Axo = 30 nm. 
Shocks form near the focal region. The simulation was done 
with about 350 points per wavelength. 
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Fig. 20. Focal pressure for an oscillation Axo = 30 nm. The 

simulation was done with about 350 points per wavelength. 



pressure does not depend on the value of the maximal 
pressure. Thus it is still possible to calculate a minimum 
pressure value and compare with experiments. 

Fig. |2l| shows the amplitude of the transducer oscilla- 
tion which is necessary to reach a given minimal pressure 
at the center, starting from a static pressure Pst- Actually, 
the variable we plot is rather Pgt as a function of pst-^ajQ, 
in order to see the departure from the linear theory 



^min — ^st + La^pPst^XQ. 



(28) 



Our interest in such a curve came from a first version of 
the experiment, for which it was not possible to measure 
the focal pressure. Still, we were able to measure the oscil- 
lation amplitude necessary to obtain cavitation for various 
initial static pressures, i.e. we could plot a curve such as 



those of Fig. Ell, for the special case Pni 



the 



cavitation pressure. However this could be done only for 
positive static pressures, while we were interested in the 
zero amplitude value of the curve, for which static pressure 
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Fig. 21. Static pressure, as a function of the oscillation am- 
plitude for which the minimal focal pressure is equal to a given 
value Pmin, multiplied by the static density. Each curve corre- 
sponds to a different value of Pmin- 



and cavitation pressure should be the same. Thus simula- 
tions are useful to determine which kind of extrapolation 
should be used for negative static pressure values. 

It is interesting to see that the effect of non-linearities 
is to bend such curves in a concave way (Fig. ^). In 
Ref. 1^ , the sign of this curvature was used to show that a 
linear extrapolation provides an upper bound of the cav- 
itation pressure. Whether the shock formation affects the 
nucleation mechanism is an open question. If we try to plot 
the same oscillation amplitude as a function of the static 
density in the cell (Fig. instead of the static pres- 
sure (Fig. non-linearities are even more pronounced, 
as would have been expected from the equation of state (a 
concave function of a concave function is still more con- 
cave) . 

More important for the validation of our theoretical 
methods is that we succeeded later in measuring at the 
focal point the temporal signal itself. As explained below, 
this was done in a quasi-spherical geometry and very good 
agreement was found between theory and experiments. 



3 Experiments 



A hemispherical piezoelectric transducer is held against a 
clean glass plate. In a first approximation, the glass re- 
flects the sound wave so that this is equivalent to a full 
spherical geometry. The main interest of the glass plate is 
that it allowed us to measure the instantaneous density at 
the center from the reflection of light at the glass/helium 
interface. Indeed, the reflectance depends on the refrac- 
tion index of liquid helium which depends on its density 
as is well known from the Clausius-Mossoti relation. 
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Fig. 22. Same figure as Fig. except that the static pressure 
is replaced by the corresponding static density. The curvature 
of the function is enhanced by using this variable. 




Fig. 23. The experimental set-up which is immersed in liq- 
uid helium, inside the experimental cell. At the center of the 
hemispherical transducer, the amplitude of the sound oscilla- 
tion is measured from the intensity of the reflected light. The 
transmitted light is used to detect the possible nucleation of 
bubbles or crystallites. 



3.1 Experimental method 

Our cxpci-imcntal method is described in full details else- 
where [ p^[20[] . Let us only summarize it here. The trans- 
ducer radius is 8 mm and its thickness is 2 mm. It resonates 
in a thickness mode at / = 1.019 MHz. At this frequency, 
it has a minimum impedance Z = 22 Ohm. Its quality 
factor is Q = 50 ± 5 when immersed in liquid helium at 
25 bar. We usually pulse it with bursts of 6 oscillations. 
The 300 cm'^ experimental cell is full of liquid helium and 
attached to the mixing chamber of a dilution refrigerator. 
We can run the experiment between 30 mK and 1.5 K, at 
static pressures from to 25 bar. Fig. |2^ shows our optical 
setup: a brass piece holds a lens whose focal length is 21.8 
mm in liquid helium and a wedged glass plate (20 mm 
in diameter, « 2 mm thick). The transducer is pressed 
against the plate. All the space inside is filled with liquid 
He. Thanks to the lens, the radius of the laser waist is re- 
duced from 320 fim to 7 fim. This means that the spatial 
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resolution is about 14 //m, the diameter of the optical focal 
region. This is small compared to the size of the acoustic 
focal region which is set by the acoustic wavelength at 
1 MHz: from 240 fim at bar to 360 fim at 25 bar. The 
distance of the lens to the glass plate has been carefully 
adjusted to have the laser focused at the glass/helium in- 
terface, on the transducer side. This was checked from the 
parallelism of the reflected beam. The 2 degrees wedge of 
the glass plate avoids interferences with reflections on its 
front face. A 1.7 mm hole in the transducer allows the 
transmitted light to be analyzed on the other side of the 
cryostat (see Fig. p4| ). 

We use a single mode Ar+ laser. It is operated around 
10 mW where its stability is best. In front of the laser is an 
absorber which reduces the light entering the cell (Fig. p^ . 
We use the reflection on the front face of the absorber to 
monitor possible drifts of the laser power, but its stability 
is better than 0.5 % per day. For low temperature mea- 
surements, we reduced the dissipation in the cell with a 
small electromechanical shutter; it was synchronized by 
the pulse generator so that the cell was illuminated dur- 
ing 10 msec only, around the arrival time of each acoustic 
pulse on the glass plate. Two moving mirrors allow us to 
translate the laser beam vertically and horizontally. The 
last mirror is mounted on a rotating plate, so that the an- 
gle of incidence of the beam can be adjusted as well. These 
rotations correspond to translations of the optical focus in 
the focal plane of the lens. By successive operations, the 
laser spot was brought to the center of the acoustic focal 
region. The final adjustment was obtained by maximiz- 
ing the modulation of the reflected signal by the acoustic 
wave. 

The transmitted light is collected by a photomulti- 
plier and used to detect nucleation events one by one. 
The light which is reflected at the glass/helium interface 
is separated from the incident beam by means of a semi- 
transparent plate and directed towards a photodiode. We 
used either a Hamamatsu C5331-03 avalanche photodiode 
(the "APD"), for the detection of the ac-modulation of the 



reflected beam, or a Hamamatsu S1406 silicon photodiode 
(the"SPD") for the detection of the dc-component. The 
output from the photodiodcs is digitized with a LeCroy 
9344 CM oscilloscope at 1 GS/s with 8-bit resolution (6.5 
effective bits due to the clock jitter at 1 GS/s). 

The ac-component of the signal is related to the mod- 
ulation of the helium density by the acoustic wave, and it 
is at most a few percent of the dc part which is related 
to the static density. The acoustic transmission into the 
glass being very small, we have neglected its effect on the 
reflected light. In order to achieve a 1% accuracy on the 
acoustic wave amplitude, we reduced the noise by averag- 
ing on 10000 sound bursts with a repetition rate of 1 to 
10 Hz. There are two main sources of noise. The first one 
is photon noise in the reflected light. For a 5 /iW power, 
this quantum noise is typically 10 nW, much more than 
the resolution we need. The other noise source is the os- 
cilloscope jitter and it has a comparable amplitude. After 
averaging, the signal is well enough extracted from the 
noise as shown on Fig. k1. 



3.2 Calibration 

The intensity of the reflected light is proportional to the 
normal reflectance R at the glass/helium interface: 



R 



(29) 



where Ug = 1.5205 is the refractive index of glass for 514.5 
nm green light. As for the index n of helium it is given by 
the Clausius Mossoti relation: 



-1 



3AI 



(30) 



where M — 4.0026 g, p is the helium density, and aM = 
0.1245 cm"^mol"^ is the molar polarizability for the same 
green light. Note that a a/ increases slightly as a func- 
tion of frequency from its zero frequency value a mo = 
0.1233 cm-^mol"^, as explained by successive authors [ pl[ . 

We proceed as follows. We flrst measure the static pres- 
sure in the cell. Knowing the equation of state P{p), we 
obtain the static density js). From the static density, we 
calculate the normal reflectance in the absence of modula- 
tion by the wave. This is our reference. We then measure 
the ratio of the ac- to the dc-component of the reflected 
light, and obtain the amplitude of the density modulation 
in the acoustic wave. Unfortunately, we cannot do this 
with a single diode, so that we have to calibrate the ra- 
tio of the respective gains of our two photodiodes. This 
is achieved in the overlap of their bandwidths, with the 
help of an acousto-optic modulator operated at 200 kHz. 
We also have to check the linearity of the APD. Its gain 
is found constant up to incoming powers of 6 pW where 
it starts decreasing. Most of the time, we use the APD in 
its linear regime; otherwise, a small correction is applied. 

We flnally have the following sources of uncertainties: 
the static density can be known within 2 to 4 10~^ kgm~^. 



ising of a spherical acoustic wave 
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32.5 33 
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Fig. 25. Two recordings of sound wave amplitudes respec- 
tively corresponding to excitation voltages 9.05 and 20.4 V on 
the transducer. The static density is p = 0.15851 g/cm^ cor- 
responding to a static pressure Pstat = 9.80 bars (horizontal 
line). The asymmetry of the oscillations is well reproduced by 
the numerical calculations performed for Axq = 3 and 6 nm. 
The numerical focal density (solid lines) is obtained with a 
spatial average weighted by a Gaussian with waist 7/im. Sim- 
ulations are performed with 350 mesh points per wavelength. 

About the same uncertainty comes from the determina- 
tion of the base hne of the APD signal. Uncertainties in 
the gain ratio and in the APD measurement lead to a 1% 
uncertainty in the wave amplitude. We could finally check 
our calibration by studying heterogeneous nucleation of 
bubbles on the glass plate. We observe various nucleation 
mechanisms. One of them occurs at saturated vapor pres- 
sure (Psv = Obar in the low temperature limit) where 
the liquid density is 145.13 kgm~^. In a series of measure- 
ment at a static pressure Pst — 4.30 bar, we found cavita- 
tion at 145.15 kgm^^; in another series of measurements 
at Pst = 2.95 bar, we found cavitation at 145.12 kgm~^. 
This illustrates the final uncertainty in our measurements. 



3.3 Comparison with calculations 



Fig. 25 shows two recordings obtained at O.IK with re- 
spective excitation amplitudes of 9.05 and 20.4 V on the 
transducer. In the cell, the static pressure is 9.80 bar, cor- 
responding to a static density p — 158.51 kgm^^. The 
density oscillation is found asymmetric at large ampli- 
tude: negative swings are broader with a smaller ampli- 
tude than positive swings. Moreover, the negative swings 
are not symmetric in time. We have chosen this recording 
at intermediate pressure and moderate amplitude because 
the signal shape is not modified by any nucleation of crys- 
tals or bubbles. 

We can compare the experimental recordi ngs with the 
numerical calculations described in sections 2.3 and 2.4. 



Simulations are performed with the same static pressure as 
in the experiment. Then there is only one free parameter 
in the simulation (the oscillation amplitude Axq) in order 
to adjust both the amplitude and the shape of the sig- 
nal. The adjustment with the experimental signal is made 




10 20 30 

Experimental Ax^ (nm) 

Fig. 26. Oscillation amplitude Axo for which the simulation 
fits the experimental signal, as a function of the experimental 
applied voltage. These results are obtained for different static 
pressures, i.e. 110 mbars (★), 800 mbars (o), 4.3 (□), 5.1 {(}), 
9.8 (A), 15.4 (V), 22 (l>), and 25.3 (+) bars. We also indicate 
as a second x-axis the estimate for the oscillation amplitude 
given by (^) which is obtained from independent physical ar- 
guments. 



only on the central oscillation of the latter. Indeed, in the 
experiment, the transducer is excited with an electrical 
burst of six oscillations. Since the transducer has a finite 
quality factor Q sa 50, the amplitude of the sound wave in- 
creases during six periods and slowly decreases afterwards. 
The numerical result is used in the steady regime, after 
the initial transient. As can be seen, we find a very good 
agreement for numerical oscillation amplitudes Axq = 3 
and 6 nm: both the asymmetry with respect to the hori- 
zontal axis and the asymmetry in time are well reproduced 
by the calculation. 

We make similar adjustments with several recordings, 
for different oscillation amplitudes and different static pres- 
sures. Our results are summarized on Fig. For each 
adjustment, the numerical amplitude obtained by fitting 
the central oscillation is associated to the experimental 
voltage applied to the transducer. On Fig. Ep, we also in- 



dicate as a second x-axis the estimate given by ( 38 ) for the 
oscillation amplitude Axq . As this estimate is obtained in- 
dependently from the simulations, as this will be detailed 
below, we refer to it as the experimental oscillation am- 
plitude in the figure. If the estimate and the numerical 
adjustment were in perfect ag reement, one would expect 
a slope equal to one in Fig. g6[ Actually, we find that both 
methods give different results, and this will be discussed 
below. The adjustment between experimental signals and 
numerical simulations yields the following calibration : 



V 



0.30 ± 0.02 nmV" 



(31) 
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It is interesting to compare this value with an estimate 
from the measurement of the electrical characteristics of 
the ceramic. Let us summarize the derivation of this esti- 
mate, which is given in Ref. [§. 

Indeed, the quality factor is simply related to the ra- 
tio of the acoustic energy _Eac which is stored during one 
period to the average dissipated power V"^/(2Z): 



(32) 



and the acoustic energy can be evaluated as follows. Let us 
call R the mean radius of the transducer and 2e its thick- 
ness. For a resonance in a thickness mode, one can assume 
that the sound wave inside the transducer is a spherical 
wave with an amplitude proportional to 1/r sin[(fc(r — 
R)] sui(u;t) where the wavevector k = 7r/2e and to = 27r/. 
For a spherical wave, the elastic displacement inside the 
transducer is 



R 

Uo — sin[fc(r — R)] s'm{ujt). 



(33) 



Since i?ac is twice the average kinetic energy in the 
transducer and the local velocity is simply the time 
derivative of u, one can integrate over the thickness and 
write: 



TrptU! UqR e 



(34) 



where pt is the transducer density. After expressing i?ac in 
terms of the maximum displacement Axq = uqR/ [R — e) 
of the inner surface, we obtain 



M 



(35) 



where M — 7.510^'^ kg is the mass of our transducer. Note 
that the last factor was forgotten in Ref. ||]. Finally, we 
can express the displacement Axq as a function of the 
applied voltage: 



Axo = V 



2Q 



1/2 



R-e 
R 



1 



3R^ 



The above equation leads to Axq/V — 1.7nmV^^. 
However, this value would correspond to an excitation 
with long bursts, when the stored energy saturates. Since 
we excite it with bursts of six oscillations only, and since 
the phase is such that they start with positive swings, the 
maximum pressure is reached (6-1-1/4) periods after time 
zero, and the displacement to be considered in our case is 

Axo - ^a^ooo [1 - exp(-12^/Q)] [exp(-7r/20)] (37) 

We find Axq ~ 0.513 Axq^ with Axq^ given by Eq. 
and our final prediction is 



^ =0.87nmV-i 



(38) 



The above value has the right order of magnitude but 
it is three times more than given by the fit of our numeri- 
cal calculations (Eq.|l]) . There are several assumptions in 
the above analysis which can be claimed as responsible for 
this discrepancy. We list them starting with those that we 
expect to be the most relevant: 

a - the resonance in a thickness mode may be coupled to 
flexion modes, in which case the efficiency of the trans- 
ducer can easily be reduced. 

b - there is a small hole in the center of the transducer 

which allows the transmitted light to be analyzed. 

c - the sound wave in the transducer cannot be strictly 

spherical, since there must be edge effects near its free 

equator. 

d - the reflexion by the glass plate is not perfect so that 
our closed hemispherical geometry is not strictly equiva- 
lent to a full spherical geometry. Once more this should 
reduce the efficiency of the transducer, 
e - some of the emitted energy is lost in the various pieces 
which hold it in the cell. 

f - there is also some uncertainty of order 10 % in the 
measurement of Z and Q. 

As a result, we consider the value 0.3 nm V^^ as a very 
useful calibration of the efficiency of our transducer, in 
qualitative agreement with a simple estimate. 



4 Conclusions and perspectives 

In this article, we have presented analytic and numerical 
calculations of the focusing of a spherical acoustic wave. 
We have shown that shocks are generated in this geometry 
and we have obtained an analytic estimate for the shock 
length based on the characteristics method. Then, in order 
to perform full numerical simulations of the focussing pro- 
cess, we have used a WENO scheme to treat shocks. We 
then showed that our method is validated by a compari- 
son with experimental measurements in a quasi-spherical 
geometry. We have measured a wave distortion which is 
well reproduced by our calculation and the analysis of its 
dependence on the excitation amplitude has led us to a 
very useful calibration of the efficiency of our transducers. 

We consider this work as a first step only, and we plan 
to extend it to a hemispherical geometry for two impor- 
tant reasons. Indeed, in order to study the homogeneous 
nucleation of bubbles in stretched fluids or that of crys- 
tals in pressurized fluids, we need to eliminate the effect of 
walls. This is achieved by using hemispherical transducers 
which focus acoustic waves away from any walls [0,^. In 
such experiments, since we have no probe in the acous- 
tic focal region where nucleation takes place, there is a 
difficult problem of calibration of the sound amplitude, 
for which any reliable calculations would be very useful. 
Of course, the calculation in a hemispherical geometry is 
much more difficult because it is two-dimensional (it de- 
pends on both the radial distance and the polar angle). 
Now that the method is known for the treatment of shocks, 
the 2-D calculations should be tried. Furthermore, we have 
estimated the amplitude of non-linear effects in the hemi- 
spherical geometry Q , and found them much smaller than 
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in the spherical geometry, though lower pressures seem 
to be reached. This is interesting in itself and should be 
tested numerically. One physical explanation could be that 
the local condition at the center is different: by symme- 
try, the spherical geometry imposes that the center is a 
node for the fluid velocity. In the hemispherical geometry, 
there is no reason why it should be so. On the contrary, 
the sound wave could even create a flow at the center with 
non vanishing averaged value. This phenomenon is known 
in the literature as acoustic streaming. This symmetry 
difference might lead to a different amplitude for the non- 
linear effects. It would be very interesting to study this 
phenomenon numerically. 

Another direction of research would deal with the in- 
teraction between shocks and nuclcation. Until now, all 
theories predicting the nucleation threshold completely 
ignore the presence of very steep gradients. This is not 
necessarily justified. 
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Appendix : theoretical prediction of the shock 
length by the characteristics method 



Appendix A: Characteristic equations 

The system of Euler equations (||) is of the form 

dtVi -f AijdrVj = hi, (39) 
with V = {p, u) and b = {—6pu/r, 0). The matrix 




has two eigenvalues p+ = u+Cg and /i_ — u — Cg associated 
with the left eigenvectors 



1. 



1 

1 



(40) 



If we apply \k on the left of equation (B9h, we obtain 



dt 



= 0, 



where 



dt 



= dt+ Hkdr 



(41) 



(42) 



Thus the derivative d/dt is taken along a curve r — x{t) 
with slope dx/dt — fik everywhere {pk itself being a func- 
tion of r and t via v). By definition, this curve is called 
a k-th characteristic and is denoted by C''^-'. There is a 
whole family of C'-'^' characteristics, covering the whole 
space, each curve being determined, for example, by the 
initial conditions. 

The leftmost term of equation (^) can be integrated 
and the equation becomes 

^ [C,{p -\np)± u] + e^^^^P^ ^ 0, (43) 
dt r 

where again the time-derivative is taken along a charac- 
teristic curve. 

We call Rieman invariants the quantities 



I± = Co{p -~\np) ±u 
appearing in eq. m2 



(44) 
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Appendix B: Lower bound for the shock length 



Considering the geometry and notations of Sec. |2.2.2 , 
we calculate at what time the first characteristic Cq 
emitted by the piston at t = is cut by another charac- 
teristic cl~\ This is the signature of a shock forming. The 
calculation assumes that it is the first shock ever formed 
in the cell. 

As Cq"-* is first cut by characteristics C^"-' emitted at 

early times, and almost parallel to Cg \ we perform the 
following change of variable. Instead of r, the location of 
any point will be given by its distance 77 to a point moving 

on the characteristic Cq \ taken at the same time : 
1] = Lq - Cstt - r. 



(45) 



The characteristic equations (^) have to be written 
in the new coordinates {ri,t). Besides, we would like to 
eliminate u and p from the equations so that the only 
remaining unknowns would be the Rieman invariants /+ 
and Following 0, we have drj/dt — — (m ± Cg + Cst) 
and u — (/_). — IS) 12. But unlike Ref. 0, the density p 
-and thus the sound velocity- cannot be expressed in a 
simple way in terms of /+ and /_ . Indeed, we would need 
to invert the relation 



/++/_ =2Co(p-lnp) = 2/(p). 



(46) 



At this stage, we will only assume that the inversion can 
be performed and write 



-1 (U^i- 



p^x- 

The characteristic equations now read 



L - Cstt - V 



Cst > drfl^ 



Co 



Then if we write equations (^8[|49|) at lowest order, all the 
terms in the second equation vanish. In the first equation, 
both the first and last terms disappear. The remaining 
term leads to 



0. 



(51) 



This is not very surprising, as corresponds to the char- 
acteristics moving from the fluid at rest into the perturbed 
region. Taking (MS) to the next order, we obtain an equa- 



tion for / 



(1) 



dl 



(1) 



dt t~Lo/cst 



r(i) 



where 



1 fipst-l 



2 \Pst-l 



0, (52) 



(53) 



This equation can be solved using a change of variables 

ilds 

/i^)(0) 



/_ = l//i^^ The solution yields 



/«(r) 



with 



(1-r) l-/Ci^/i^^(0) In(l-T) 



, (54) 



(55) 



We have now to determine the initial value /i^'' (0). It refers 
to small t, rather than t exactly equal to zero. 

To determine its value, we calculate for small time t 
the variation of /_ between point A = (r, t) = {LQ — Cstt, t) 

(47) which sits on the first characteristic Cq •* , and point B = 
{r,t) = {rp{t),t) where rp{t) is the location of the piston 
at time t. 

At A, the fiuid is at rest. We have /_ = /i"' and rj ^ 0. 
At B, for t small, r ~ Lq and thus rj ~ —Cstt- Besides, 
/_ — Co{p — hi p) — u where u is equal to the velocity of 
the piston Vp{t) = —Avq sm{ujt) ~ —AvQCut. If we expand 
p{rp{t),t) = pst + at (where a is unknown), then replacing 
into /_ and expanding in t yields 



I - 1 



= 
(48) 



L - Cstt - ri 



Co 



/_ = Coipst - hipst) + Coat - Co — t + AvoLut 

Pst 

^ /(") + (^^a + Avoio^ t. 



= 

(49) Comparing this relation with the expansion 



We search for a solution under the form 

00 

/4^,i) = E^+'^W^"' 



m=0 



m=0 



(50) 



We choose the lowest order terms /i"'' and /i^-* equal to 



their value in the fluid at rest / 



(0) _ r(0) 



Co(pst-lnpst). 



77/ 



(1) 



AO) 



CstlL\ 



gives 



(56) 



^ _^ _ (57) 

Pst Cst 

On the other hand, an expansion of ( p7| ) in powers of rj 
gives 



at = PstPst = — V^r- t 
Pst \ 2 / 



(58) 
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The elimination of a between the above two equations 
yields 



The time at which /i^^ becomes infinite (see equation 
|5^ ) gives an upper bound tshock for shock formation. It 
is only an upper bound because some other terms of the 
expansion ([50|) in rj may explode before Wc find 



Cst 



(59) 



t 



shock 



< 



Lo 

Cst 



1 — exp 



Pst 



1 



2LoujAvo 



< 



Lo 

Cst' 



(60) 

As the corresponding shock distance rshock is measured 
from the center of the sphere, a lower bound for rshock is 



Lq > r-shock > Lq exp 



Pst - 1 



2Lau;Avo pst 



>0. (61) 



All the above calculations are valid for small rj, i.e. 
only for characteristics not too far from Cq These are 
the characteristics emitted by the initial motion of the 
piston. 



